function u=couple60(N,SNRx,SNRy,delay_x,delay_y,filter_flag)
%
% Input: SNRx, SNRy - for each channel
%        delay_x,delay_y - 60 hertz delay in samples. 
%
% Output: u - 2 x N time series matrix
%

%
% Example generation
%

%
% signal to noise ratios
%SNRx=40;
%SNRy=40;
% Coupling parameter
beta=2;
%
Nx=3000; % burn-in length to encompass comb filter
%
%
%N=1024; % number of data points
e=rand(2,N+Nx);
x=zeros(N+2+Nx,1);
y=x;
% Dynamical coupling
for k=3:N+2+Nx;
    x(k)=1.5806*x(k-1)-.7225*x(k-2)+e(1,k-2);
    y(k)=.5*y(k-1)+e(2,k-2)+ beta*x(k-1);
end
x=x(3:N+2+Nx);
y=y(3:N+2+Nx);
Nt=Nx+2+N-3+1;
% square root powers
sx=std(x);
sy=std(y);
%
Ax=sx*10^(-SNRx/20);
Ay=sx*10^(-SNRx/20);
%
%
dx=zeros(delay_x+1);
dx(delay_x+1)=1;
%
%
dy=zeros(delay_y+1);
dy(delay_y+1)=1;
%

xu=filter(dx,1,Ax*cos(2*pi*60*(0:Nt-1)/250));
yu=filter(dy,1,Ay*cos(2*pi*60*(0:Nt-1)/250));
%xu=zeros(size(x))';
%yu=zeros(size(y))';
% final signal
xt=x+xu';
yt=y+yu';
%xt=xu';
%yt=yu';
%
xt=xt-mean(xt);
yt=yt-mean(yt);
if filter_flag==1
% notch filter parameters for BW=.5 Hz - 250 Hz sampling for 60 Hz.
 num=[ 0.9978   -0.1253    0.9978];
% roots(num)
 den=[1.0000   -0.1253    0.9956];
% roots(den)
% %
% figure
% plot(roots(den),'b+')
% hold
% plot(roots(num),'ro')
% axis([-1 1 -1 1])
% axis('square')
% hold
%pause
%
xf=filter(num,den,xt);
yf=filter(num,den,yt);
elseif filter_flag==2
%
% regressive filtering
%
[xx,xf,aax,bbx]=ad_sinus(xt,60,250);
[yy,yf,aay,bby]=ad_sinus(yt,60,250);
% aax
% bbx
% aay
% bby
%pause
%
else
xf=xt;
yf=yt;
end
%
% figure
% plot(x)
% hold
% plot(y,'r')
% hold
% figure
% plot(xu)
% hold
% plot(yu,'r')
% hold
% figure
% plot(xt)
% hold
% plot(yt,'r')
% hold
% title('t')
% figure
% plot(xf)
% hold
% plot(yf,'r')
% hold
% title('f')
% %
% %keyboard
%
xf=xf(Nx+1:Nx+N);
yf=yf(Nx+1:Nx+N);
%
u=[xf yf];
